{
 "metadata": {
  "signature": "sha256:dc92ea161c1f427380ed83be00655f8de47c597737a2b73085c20ea33aced2e9"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Theoretical ecology: Hastings and Powell\n",
      "========================================\n",
      "\n",
      "Overview\n",
      "--------\n",
      "\n",
      "A simple script that recreates the min/max bifurcation diagrams from\n",
      "Hastings and Powell 1991.\n",
      "\n",
      "Library Functions\n",
      "-----------------\n",
      "\n",
      "Two useful functions are defined in the module bif.py."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "\n",
      "def window(data, size):\n",
      "    \"\"\"A generator that returns the moving window of length\n",
      "    `size` over the `data`\n",
      "\n",
      "    \"\"\"\n",
      "    for start in range(len(data) - (size - 1)):\n",
      "        yield data[start:(start + size)]\n",
      "\n",
      "\n",
      "def min_max(data, tol=1e-14):\n",
      "    \"\"\"Return a list of the local min/max found\n",
      "    in a `data` series, given the relative tolerance `tol`\n",
      "\n",
      "    \"\"\"\n",
      "    maxes = []\n",
      "    mins = []\n",
      "    for first, second, third in window(data, size=3):\n",
      "        if first < second and third < second:\n",
      "            maxes.append(second)\n",
      "        elif first > second and third > second:\n",
      "            mins.append(second)\n",
      "        elif abs(first - second) < tol and abs(second - third) < tol:\n",
      "            # an equilibrium is both the maximum and minimum\n",
      "            maxes.append(second)\n",
      "            mins.append(second)\n",
      "\n",
      "    return {'max': numpy.asarray(maxes),\n",
      "            'min': numpy.asarray(mins)}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Model\n",
      "---------\n",
      "\n",
      "For speed the model is defined in a fortran file and compiled into a\n",
      "library for use from python. Using this method gives a 100 fold increase\n",
      "in speed. The file uses Fortran 90, which makes using f2py especially\n",
      "easy. The file is named hastings.f90."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "```\n",
      "module model\n",
      "    implicit none\n",
      "\n",
      "    real(8), save :: a1, a2, b1, b2, d1, d2\n",
      "\n",
      "contains \n",
      "\n",
      "    subroutine fweb(y, t, yprime)\n",
      "        real(8), dimension(3), intent(in) :: y\n",
      "        real(8), intent(in) :: t\n",
      "        real(8), dimension(3), intent(out) :: yprime\n",
      "\n",
      "        yprime(1) = y(1)*(1.0d0 - y(1)) - a1*y(1)*y(2)/(1.0d0 + b1*y(1)) \n",
      "        yprime(2) = a1*y(1)*y(2)/(1.0d0 + b1*y(1)) - a2*y(2)*y(3)/(1.0d0 + b2*y(2)) - d1*y(2)\n",
      "        yprime(3) = a2*y(2)*y(3)/(1.0d0 + b2*y(2)) - d2*y(3)\n",
      "    end subroutine fweb\n",
      "\n",
      "end module model\n",
      "```"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Which is compiled (using the gfortran compiler) with the command:\n",
      "``f2py -c -m hastings hastings.f90 --fcompiler=gnu95``"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "from scipy.integrate import odeint\n",
      "import bif\n",
      "\n",
      "import hastings\n",
      "\n",
      "# setup the food web parameters\n",
      "hastings.model.a1 = 5.0 \n",
      "hastings.model.a2 = 0.1 \n",
      "hastings.model.b2 = 2.0 \n",
      "hastings.model.d1 = 0.4 \n",
      "hastings.model.d2 = 0.01\n",
      "\n",
      "# setup the ode solver parameters\n",
      "t = numpy.arange(10000)\n",
      "y0 = [0.8, 0.2, 10.0]\n",
      "\n",
      "def print_max(data, maxfile):\n",
      "    for a_max in data['max']:\n",
      "        print >> maxfile, hastings.model.b1, a_max\n",
      "\n",
      "x_maxfile = open('x_maxfile.dat', 'w')\n",
      "y_maxfile = open('y_maxfile.dat', 'w')\n",
      "z_maxfile = open('z_maxfile.dat', 'w')\n",
      "for i, hastings.model.b1 in enumerate(numpy.linspace(2.0, 6.2, 420)):\n",
      "    print i, hastings.model.b1\n",
      "    y = odeint(hastings.model.fweb, y0, t)\n",
      "\n",
      "    # use the last 'stationary' solution as an intial guess for the\n",
      "    # next run. This both speeds up the computations, as well as helps\n",
      "    # make sure that solver doesn't need to do too much work.\n",
      "    y0 = y[-1, :]\n",
      "\n",
      "    x_minmax = bif.min_max(y[5000:, 0]) \n",
      "    y_minmax = bif.min_max(y[5000:, 1]) \n",
      "    z_minmax = bif.min_max(y[5000:, 2]) \n",
      "\n",
      "    print_max(x_minmax, x_maxfile)\n",
      "    print_max(y_minmax, y_maxfile)\n",
      "    print_max(z_minmax, z_maxfile)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}